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ABSTRACT 

The development of benchmark examples for quasi-static delamination propagation 
prediction is presented. The example is based on a finite element model of the Mixed- 
Mode Bending (MMB) specimen for 50% mode II. The benchmarking is 
demonstrated for Abaqus/Standard®, however, the example is independent of the 
analysis software used and allows the assessment of the automated delamination 
propagation prediction capability in commercial finite element codes based on the 
virtual crack closure technique (VCCT). First, a quasi-static benchmark example was 
created for the specimen. Second, starting from an initially straight front, the 
delamination was allowed to propagate under quasi-static loading. Third, the load- 
displacement as well as delamination length versus applied load/displacement 
relationships from a propagation analysis and the benchmark results were compared, 
and good agreement could be achieved by selecting the appropriate input parameters. 
The benchmarking procedure proved valuable by highlighting the issues associated 
with choosing the input parameters of the particular implementation. Overall, the 
results are encouraging, but further assessment for mixed-mode delamination fatigue 
onset and growth is required. 


INTRODUCTION 

Over the past two decades, the use of fracture mechanics has become common 
practice to characterize the onset and growth of delaminations. In order to predict 
delamination onset or growth, the calculated strain energy release rate components are 
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compared to interlaminar fracture toughness properties measured over a range from 
pure mode I loading to pure mode II loading. 

The virtual crack closure technique (VCCT) is widely used for computing energy 
release rates based on results from continuum (2D) and solid (3D) finite element (FE) 
analyses and to supply the mode separation required when using the mixed-mode 
fracture criterion [1,2]. The virtual crack closure technique was recently implemented 
into several commercial finite element codes such as Abaqus/Standard' MD 
Nastran™ 1 2 , Marc™ 2 and ANSYS 3 . As new methods for analyzing composite 
delamination are incorporated into finite element codes, the need for comparison and 
benchmarking becomes important since each code requires specific input parameters 
unique to its implementation. These parameters are unique to the numerical approach 
chosen and do not reflect real physical differences in delamination behavior. 

An approach for assessing the mode I, mode II and mixed-mode I and II, 
delamination propagation capabilities in commercial finite element codes under static 
loading was recently presented and demonstrated for VCCT for ABAQUS R [3,4] as 

TM TM _ ~~ 

well as MD Nastran and Marc [5], First, benchmark results were created manually 
for finite element models of the mode I Double Cantilever Beam (DCB), the mode II 
End Notched Flexure (ENF) and the mixed-mode I/II Single Leg Bending (SLB) 
specimen. Second, starting from an initially straight front, the delamination was 
allowed to propagate using the automated procedure implemented in the finite element 
software. The approach was then extended to allow the assessment of the delamination 
fatigue growth prediction capabilities in commercial finite element codes [4,6]. 

The objective of the present study was to create additional benchmark examples to 
simulate mixed-mode I/II conditions based on the Mixed-Mode Bending (MMB) 
specimen shown in Figure 1. Static benchmark results were created for three mixed- 
mode ratios (G///Gt= 0.2, 0.5, and 0.8) based on the approach developed earlier [3]. 
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Figure 1. Mixed-Mode Bending (MMB) Specimen. 


1 ABAQUS® is a product of Dassault Systemes Simulia Coip. (DSS), Providence, RI, USA 

2 MD Nastran™ and Marc™ are manufactured by MSC. Software Corp., Santa Ana, CA, USA. 
NASTRAN® is a registered trademark of NASA. 

3 ANSYS 8 is a product of ANSYS, Inc., Canonsburg, PA, USA 
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To create the benchmark results, two-dimensional finite element models were 
used for simulating the MMB specimens with different delamination lengths ao. For 
each delamination length modeled, the load, P, and the displacement, w, were 
monitored. The total strain energy release rate, GV, and mixed-mode ratio, Gi/Gt, were 
calculated for a fixed applied displacement. The delamination was assumed to 
propagate when the computed energy release rate, Gt, reaches the mixed-mode 
fracture toughness G c . Thus, critical loads and critical displacements for delamination 
propagation were calculated for each delamination length modeled. From these critical 
load/displacement results, benchmark solutions were created. It is assumed that the 
load/displacement relationship computed during automatic propagation should closely 
match the benchmark cases. 

After creating the benchmark cases, the approach was demonstrated for the 
commercial finite element code Abaqus/Standard® . Starting from an initially straight 
front, the delamination was allowed to propagate under quasi-static loading based on 
the algorithms implemented into the software. Input control parameters unique to 
Abaqus/Standard R were varied to study the effect on the computed delamination 
propagation. These parameters included the release tolerance and stabilization factors 
as discussed later. The benchmark enabled the selection of the appropriate input 
parameters that yielded good agreement between the results obtained from the 
propagation analysis and the benchmark results. Once the parameters have been 
identified, they may then be used with confidence to model delamination growth for 
more complex configurations. 

In this paper, only the development of one benchmark case ( Gi/Gt =0.5) for the 
assessment of the quasi-static delamination propagation is presented. Examples of 
automated propagation analyses are shown, and the selection of the required code 
specific input parameters are discussed. More detailed infonnation about the other two 
benchmark cases for mixed-mode ratios Gi/Gt =0.2 and 0.8 may be found in 
reference 7. 


METHODOLOGY BASED ON FRACTURE MECHANICS 

For the current numerical investigation, the Mixed-Mode Bending (MMB) 
specimen, as shown in Figure 1, was chosen since it is simple, and exhibits the mixed- 
mode I/II opening fracture mode over a wide range of mixed-mode ratios Gi/Gt. For 
this paper, only results from finite element models for one mixed-mode ratios ( Gi/Gt 
=0.5) are discussed. Fracture mechanics concepts [8], were applied to the MMB 
specimen to create the benchmark example. For the current study, MMB specimens 
made of IM7/8552 graphite/epoxy with a unidirectional layup, [0] 2 4, were modeled. 
The material, layup, overall specimen dimensions including initial crack length, ao, 
were identical to specimens used in related experimental studies [9]. The material 
properties are given in Table I. 


TABLE I. MATERIAL PROPERTIES [7], 


E\\ = 161 GPa 

£'22 = 11.38 GPa 

£33 = 11.38 GPa 

vi 2 = 0.32 

vi 3 = 0.32 

V 23 = 0.45 

Gi 2 = 5.2 GPa 

G 13 = 5.2 GPa 

G23 = 3.9 GPa 
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A quasi-static mixed-mode fracture criterion is determined by plotting the 
interlaminar fracture toughness, G c , versus the mixed-mode ratio, Gu/Gt as shown in 
Figure 2. The fracture criterion is generated experimentally using pure Mode I 
(G//Gy=0) Double Cantilever Beam (DCB) tests, pure Mode II (G///Gy=l) End- 
Notched Flexure (ENF) tests [9], and Mixed Mode Bending (MMB) tests of varying 
ratios of Gy and Gu . For the material used in this study, the experimental data (open 
symbols) and mean values (filled symbols) are shown in Figure 2. A 2D fracture 
criterion was suggested by Benzeggah and Kenane [10] using a simple mathematical 
relationship between G c and Gu/Gt 


(Gii\ V 

G c — G lc + (G//c - G lc ) (1) 

In this expression, typically called B-K criterion, G* and G// c are the experimentally 
determined fracture toughness data for mode I and II as shown in Figure 2. The 
exponent 77 was detennined by a curve fit using the Levenberg-Marquardt algorithm in 
the KaleidaGraph™ graphing and data analysis software [11]. The parameters G/ c , G// c 
and 77 are required input to perfonn a VCCT analysis in Abaqus/Standard Rj , as 
discussed in detail in reference 7. 

During an automated propagation analysis, the total strain energy release rate, Gy, 
and the mixed-mode ratio Gu/Gt are computed using VCCT. The failure index, Gy/G c , 
is calculated by correlating the computed total energy release rate, Gy, with the mixed- 
mode fracture toughness, G c , of the graphite/epoxy material. As shown in equation (1), 
the mixed-mode fracture toughness, G c , is a function of the mixed-mode ratio Gu/Gt 
(see also Figure 2). It is assumed that the delamination propagates when the failure 
index, Gy/G c , reaches unity. 



mixed mode ratio GfG T 


Figure 2. Mixed mode fracture criterion for IM7/8852. 
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FINITE ELEMENT MODELING 


Model description 

An example of a two-dimensional finite element model of a Mixed-Mode Bending 
(MMB) specimen with boundary conditions is shown in Figure 3 for a mixed-mode 
ratio Gu/Gt =0.5. Based on previous experience [3,4,6], the specimen was modeled 
with solid plane strain elements (CPE4I) in Abaqus/Standard R 6.10. Along the length, 
all models were divided into different sections with different mesh refinement as 
shown in Figure 3a. The MMB specimen was modeled with six elements through the 
specimen thickness (2h) as shown in the detail of Figure 3b. The resulting element 
length at the delamination tip was Aa=0.5 mm. The load apparatus was modeled 
explicitly using rigid beam elements (R2D2) as shown in Figure 3a. Multi-point 
constraints were used to connect the rigid elements with the planar model of the 
specimen and enforce the appropriate boundary conditions as shown in Figures 3b and 
c. 

The plane of delamination was modeled as a discrete discontinuity in the center of 
the specimen. For the analysis in Abaqus/Standard* 6.10, the models were created as 
separate meshes for the upper and lower part of the specimens with identical nodal 
point coordinates in the plane of delamination [12]. Two surfaces (top and bottom 
surface) were defined to identify the contact area in the plane of delamination as 
shown in Figures 3b and c. Additionally, a node set was created to define the intact 
(bonded nodes) region. 

An example of a three-dimensional finite element models of the MMB specimen 
is shown in Figure 4. Along the length, all models were divided into different sections 
with different mesh refinement. A refined mesh was used in the center of the MMB 
specimen as shown in the detail of Figure 4b. Across the width, a uniform mesh (25 
elements) was used to avoid potential problems at the transition between a coarse and 
finer mesh [3-6]. Through the specimen thickness (2h), six elements were used as 
shown in the detail of Figure 4b. The resulting element length at the delamination tip 
was Aa=0.5 mm. The specimen was modeled with solid brick elements (C3D8I), 
which had yielded excellent results in previous studies [3,4,6]. The load apparatus was 
modeled explicitly using rigid plate elements (R3D4) as shown in Figure 4a. As before 
for the two-dimensional model, multi-point constraints were used to connect the rigid 
elements with the solid model of the specimen and enforce the appropriate boundary 
conditions. 

Static delamination propagation analysis 

For the automated delamination propagation analysis, the VCCT implementation 
in Abaqus/Standard* 6.10 was used. The plane of delamination in three-dimensional 
analyses is modeled using the existing Abaqus/Standard R crack propagation capability 
based on the contact pair capability [12]. Additional element definitions are not 
required, and the underlying finite element mesh and model does not have to be 
modified [12]. The implementation offers a crack and delamination propagation 
capability in Abaqus/Standard R . During the analysis, the energy release rate at the 
crack tip is calculated at the end of a converged increment. Once the energy release 
rate exceeds the critical strain energy release rate (including the user-specified mixed- 


5 



mode criteria as shown in Figure 2), the node at the crack tip is released in the 
following increment, which allows the crack to propagate. To avoid sudden loss of 
numerical stability when the crack tip is propagated during the analysis, the force at 
the crack tip is released gradually during succeeding increments. The release is 
perfonned in such a way that the force is brought to zero before the next node along 
the crack path begins to open [12,13]. 

In addition to the mixed-mode fracture criterion, Abaqus/Standard R requires 
additional input for the propagation analysis using VCCT. If a user specified release 
tolerance is exceeded in an increment (G-G C )IG C > release tolerance, a cutback 
operation is perfonned which reduces the time increment. In the new smaller 
increment, the strain energy release rates are recalculated and compared to the user 
specified release tolerance. The cutback reduces the degree of overshoot and improves 
the accuracy of the local solution [12]. A release tolerance of 0.2 is suggested in the 
handbook [12]. To help overcome convergence issues during the propagation analysis, 
Abaqus/Standard ® provides: 

• contact stabilization that is applied across only selected contact pairs and 
used to control the motion of two contact pairs while they approach each 
other in multi-body contact. The damping is applied when bonded contact 
pairs debond and move away from each other [12.13] 

• automatic, or static, stabilization that is applied to the motion of the entire 
model and is commonly used in models that exhibit statically unstable 
behavior such as buckling [12,13] 

• viscous regularization that is applied only to nodes on contact pairs that 
have just debonded. The viscous regularization damping causes the tangent 
stiffness matrix of the softening material to be positive for sufficiently 
small time increments. Viscous regularization damping is similar to the 
viscous regularization damping provided for cohesive elements and the 
concrete material model in Abaqus/Standard R [12,13], 

Further details about the required input parameters are discussed in reference 8. 

For automated propagation analysis, it was assumed that the computed behavior 
should closely match the benchmark results created below. For all analyses, the elastic 
constants and the input to define the fracture criterion were kept constant. Based on 
poor results in a previous investigation [3], automatic or static stabilization was not 
used in this study. Viscous regularization had yielded good results in the past [3,4], 
however, it was not used in the current study to limit the number of analyses. 
Instead, the study focused on using the default values and contact stabilization 
which had recently yielded good results [4]. This approach was used to gain 
confidence in the parameters that had been identified previously and thus detennine 
if the selection of the parameters could be considered problem independent. 
Therefore, only the following items were varied to study the effect on the automated 
delamination propagation behavior during the analysis. 

• The release tolerance ( reltol ) was varied. 

• Analyses were perfonned with and without contact stabilization. For 
analyses that included contact stabilization only, a single stabilization factor 
(cs= 1 x 1 O' 6 ) was used. 

• Two- and three-dimensional models with different types of elements were 
used. 
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(a). Mixed mode ratio GfGj—0.5, c=41 .3 mm. 
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(c). Detail of load introduction zone. 


Figure 3. Deformed 2D FE-model of a MMB specimen. 
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(b). Detail of specimen tip and crack tip zone. 


Figure 4. Deformed 3D FE-model of a MMB specimen. 


DEVELOPMENT OF THE STATIC BENCHMARK CASE 

The static benchmark case was created based on the approach developed earlier [3]. 
Two-dimensional finite element models simulating MMB specimens with 16 different 
delamination lengths ao were created (25.4 mm< ao< 70.6 mm). For each delamination 
length modeled, the load, P, and displacement, w, were monitored as shown in 


Figure 5 (colored lines) for the case of 50% mode II (Gi/Gt =0.5). Using VCCT, the 
total strain energy release rate, Gt, and the mixed-mode ratio Gi/Gt were computed at 
the end of the analysis as shown in Figure 5. The failure index Gt/G c was calculated 
by correlating the computed total energy release rate, Gt, with the mixed-mode 
fracture toughness, G c , of the graphite/epoxy material. It is assumed that the 
delamination propagates when the failure index Gt/G c reaches unity. Therefore, the 
critical load, P C nt, and critical displacement, w crih can be calculated for each 
delamination length modeled as discussed in detail in reference 7. The results were 
included in the load/displacement plots as shown in Figure 5. 

These critical load/displacement results indicated that, with increasing delamination 
length, less load is required to extend the delamination up to ao ~48 mm. To further 
extend the delamination, an increase in load is required past w= 2.0 mm. For the first 
three delamination lengths, ao, plotted in Figure 5, the values of the critical 
displacements also decreased at the same time. This means that the MMB specimen 
exhibits unstable delamination propagation under load control as well as displacement 
control in this region for Gi/Gt =0.5. The remaining critical load/displacement results 
indicated stable propagation. From these critical load/displacement results (solid grey 
circles and dashed line), two benchmark solutions can be created as shown in Figure 6. 
During the analysis, either prescribed displacements, w, or nodal point loads, P, are 
applied. For the case of prescribed displacements, w, (dashed blue line), the applied 
displacement must be held constant over several increments once the critical point 
(Pent, Wait) is reached, and the delamination front is advanced during these increments. 
Once the critical path (dashed grey line) is reached, the applied displacement is 
increased again incrementally. For the case of applied nodal point loads (dashed red 
line), the applied load must be held constant while the delamination front is advanced 
during these increments. Once the critical path (dashed grey line) is reached, the 
applied load is increased again incrementally. 



applied displacement tv, mm 

Figure 5. Calculated critical load-displacement behavior for a MMB specimen. 
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The benchmark result for prescribed displacements may also be visualized by 
plotting the prescribed displacements, w, at delamination growth onset versus the 
increase in delamination length, a* as illustrated in Figure 7. For the case of applied 
nodal point loads, the benchmark result may be visualized by plotting the applied 
loads, P, at delamination growth onset versus the increase in delamination length, a *, 
as illustrated in Figure 8. 



displacement w, mm 

Figure 6. Benchmark cases for applied load, P, and displacement, w. 


applied 
displacement 
w, mm 



Figure 7. Benchmark case for applied displacement, w. 
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Figure 8. Benchmark case for applied load, P. 


This way of presenting the results may be of advantage for large structures where 
local delamination propagation may have little effect on the global stiffness of the 
structure and may therefore not be visible in a global load/displacement plot. 
However, extracting the delamination length a from the finite element results required 
more manual, time consuming, post-processing of the results compared to the 
relatively simple and readily available output of nodal displacements and forces. The 
load/displacement, load/delamination-length or applied displacement/delamination- 
length relationship computed during automatic propagation should closely match the 
benchmark cases. 


RESULTS FROM AUTOMATED PROPAGATION ANALYSIS 
Computed delamination propagation for applied displacement 

The propagation analysis was perfonned in two steps using the models shown in 
Figures 3a and 4a starting from an initial delamination length, ao=25A mm. In the first 
step, a displacement of w= 1.2 mm was applied which nearly equaled the critical 
displacement, w cr if= 1-34 mm, determined earlier. In the second step, the applied 
displacement was increased to vv=8.0 mm. For this second step, automatic 
incrementation was used in Abaqus/Standard ' and a small increment size (0.5% of the 
total step) was chosen at the beginning of the step. To avoid termination, the analysis 
was allowed to cut back to an increment size of 10‘ 18 of the total step. The analysis 
was limited to 5000 increments. Further details are discussed in reference 7. 

Initially, analyses were performed using two-dimensional planar models (shown in 
Figure 3a) without stabilization or viscous regularization. Release tolerance values 
equal to the default value {reltol= 0.2) [11] and smaller were used since these values 
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had yielded better results in the past [3,4]. Using the default value ( reltol=0.2 ), the 
computed load and displacement exceeded the critical point before the initial load 
drop occurred and delamination propagation started (solid green line) as shown in 
Figure 9 where the computed resultant force (load P ) is plotted versus the applied 
displacement w. The computed path initially stayed above the benchmark, however, as 
the displacement continued to increase over 2.5 mm, the computed load/displacement 
path converged to the benchmark result (solid grey circles and dashed grey line). To 
reduce the observed overshoot, the release tolerance was decreased. For a release 
tolerance reltol= 0.1, the overshoot was reduced and the results improved (solid blue 
line) and shifted towards the benchmark case. For a release tolerance reltol= 0.01, the 
analysis terminated once the critical point was reached (solid red line) due to 
convergence problems. 

Based on problems identified during previous analyses, automated or static 
stabilization was not used in this study [3]. The results computed when contact 
stabilization ( cs ) was added are plotted in Figure 10. Using a stabilization factor of 
cs= 1 x 1 0 , and a release tolerance (reltol= 0.01) yielded results that were in excellent 
agreement with the benchmark (solid red line). Due to the fine mesh, only a small 
saw-tooth pattern was observed. The result confirms previous observations where 
small release tolerance values (reltol= 0.01) in conjunction with a stabilization factor of 
cs= 1 x 1 0‘ 6 had yielded excellent agreement with the benchmark results [3,4]. 

An alternative way to plot the benchmark is shown in Figure 1 1 where the applied 
displacement w is plotted versus the increase in delamination length a*. The results 
plotted in Figure 1 1 are the examples that were discussed above and were shown in 
the global load/displacement plots of Figures 9 and 10. The conclusions that can be 
drawn from this plot are identical to those discussed above. 



applied displacement w , mm 


Figure 9. Computed critical load-displacement behavior obtained from 2D planar 
models (CPE4) with differen t release tolerance settings. 
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Figure 10. Computed critical load-displacement behavior obtained from 2D planar 
models (CPE4) with added contact stabilization. 


applied 
displacement 
w, mm 



increase in delamination length a* mm 

Figure 1 1 . Computed critical displacement-propagation behavior obtained from 2D 

planar models (CPE4). 


The result obtained from a three-dimensional analysis is shown in Figure 12. 
Based on the results from two-dimensional planar models shown above, contact 
stabilization was added to the analysis to help overcome convergence issues. Since the 
computer time increases with tighter release tolerances, a tolerance value (reltol= 0.1) 
was selected. 
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applied displacement w, mm 

Figure 12. Computed critical load-displacement behavior obtained from 3D solid 

models (C3D8I). 

The result (solid blue line) computed for a release tolerance value reltol= 0. 1 and a 
small stabilization factor cs=lxl0" 6 are plotted in Figure 12. The load dropped and 
delamination propagation started upon reaching the critical point of the benchmark 
solution (solid circles and dashed grey line). The load/displacement path then ran 
parallel to the benchmark result but was shifted slightly towards lower loads. 
Deviation from the benchmark may be explained by the fact that the benchmark 
results were created using two-dimensional planar finite element models of the MMB 
specimen. Three-dimensional models, which provide a better approximation of reality, 
however, yield a failure index that is somewhat higher across the entire width than the 
results from two-dimensional planar finite element models as discussed in detail in 
reference 7. Therefore, delamination propagation is expected to start prior to the 
benchmark results obtained from two-dimensional planar finite element analysis, 
which ultimately leads to a shift of the entire results plot towards lower loads. 

Computed delamination propagation for applied quasi-static load 

The propagation analysis was perfonned in two steps using the models shown in 
Figures 3a and 4a starting from an initial delamination length, aa=25A mm. In the first 
step, a load P= 360 N was applied which equaled nearly the critical load, /f ;n 7=385 N, 
determined earlier. In the second step, the total load was increased (F^bOO N). The 
settings for incrementation were identical to the values mentioned in the section on 
applied displacements. 

The same steps discussed in the section on applied displacement were followed. 
Initially, analyses were performed using two-dimensional planar models without 
stabilization or viscous regularization. For the default value of reltol= 0.2, the load 
stopped increasing after reaching the critical point, but the analysis terminated 
immediately (solid green line) due to convergence problems as shown in Figure 13. 


14 


The release tolerance was not increased as suggested by the Abaqus/Standard" error in 
the message (.msg) file. Previous analysis had shown that, by increasing the release 
tolerance (reltol> 0.2), termination of the analysis could be avoided. However, the 
results had not been in good agreement with the benchmark [3,4]. Therefore, 
additional stabilization had to be introduced in order to obtain agreement with the 
benchmark case. 

The results computed when contact stabilization ( cs ) was added are plotted in 
Figure 14. A small stabilization factor (cs=lxl0' 6 ) was used for all cases since it had 
yielded good results in previous analyses [3, 4], Initially, the release tolerance value 
was set at the default value (reltol= 0.2) (solid green line). For this parameter 
combination, the load increased up to the critical point, and delamination propagation 
started while the load remained constant (solid green line). Also for the stable 
propagation path, the result was in good agreement with the benchmark result (solid 
grey line). Then, the release tolerance was reduced to reltol= 0.1 (solid blue line) and 
further to reltol= 0.01 (solid red line). For all cases, the results were in good agreement 
with the benchmark results over the entire load/displacement range. 

An alternative way to plot the benchmark is shown in Figure 15 where the applied 
load P is plotted versus the increase in delamination length a*. The results plotted in 
Figure 15 are the same as those that were discussed above and were shown as global 
load/displacement plots in Figure 14. The conclusions that can be drawn from the 
plots in Figure 15 are identical to those discussed in the above for Figure 14. 



displacement w, mm 


Figure 13. Computed critical load-displacement behavior obtained from 2D planar 
models (CPE4I) with different release tolerance settings. 
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applied 
load 
P, N 



displacement w, mm 

Figure 14. Computed critical load-displacement behavior obtained from 2D planar 
models (CPE4I) with added contact stabilization. 



increase in delamination length a* mm 

Figure 15. Computed critical load-propagation behavior obtained from 2D planar 

models (CPE4I). 


A result obtained from a three-dimensional model is shown in Figure 16. Based on 
the results from two-dimensional planar models shown above, contact stabilization 
was added to help overcome convergence issues. A release tolerance value 
(reltol=0A) was used which had yielded good results previously. The 
load/displacement result computed when a small stabilization factor (cs=lxl0" 6 ) was 
added is plotted in Figure 16. 
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displacement w, mm 

Figure 16. Computed critical load-displacement behavior obtained from 3D solid 

models (C3D8I). 


The load increase stopped and delamination propagation started (solid blue line) 
shortly before reaching the critical point of the benchmark solution (solid grey line). 
As mentioned earlier and discussed in detail in reference 7, deviation from the 
benchmark may be explained by the fact that the benchmark results were created using 
two-dimensional planar finite element models of the MMB specimen. 


SUMMARY AND CONCLUSIONS 

The development of a benchmark example, which allows the assessment of the 
static delamination propagation prediction capabilities was presented and 
demonstrated for Abaqus/Standard R . The example is based on finite element models 
of the mixed-mode I/II Mixed-Mode Bending (MMB) specimen for 50% mode II. The 
models are independent of the analysis software used and allow the assessment of the 
automated delamination growth prediction capabilities in commercial finite element 
codes based on the virtual crack closure technique (VCCT). 

First, a quasi-static benchmark result was created based on the approach 
developed in reference 3. Two-dimensional finite element models were used for 
simulating the MMB specimen with 50% mode II and different initial delamination 
lengths, ao. For each delamination length modeled, the load and displacements were 
monitored. The mixed-mode I/II strain energy release rate was calculated for a fixed 
applied displacement. The delamination was assumed to propagate when the total 
strain energy release rate reached the fracture toughness value. Thus, critical loads and 
critical displacements for delamination propagation were calculated for each initial 
delamination length modeled. From these critical load/displacement results, 
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benchmark solutions were created. The load/displacement relationship computed 
during automatic propagation should closely match the benchmark cases. 

After creating the benchmark cases, the approach was demonstrated for the 
commercial finite element code Abaqus/Standard® . Starting from an initially straight 
front, the delamination was allowed to propagate under quasi-static loading based on 
the algorithms implemented into the software. Input control parameters were varied to 
study the effect on the computed delamination propagation and growth. 

The results showed the following: 

• The benchmarking procedure proved valuable by highlighting the issues 
associated with choosing the input parameters of the particular 
implementation. 

• In general, good agreement between the results obtained from the 
automated propagation analysis and the benchmark results could be 
achieved by selecting input parameters that had previously been 
determined during analyses of mode I Double Cantilever Beam and 
mode II End Notched Flexure specimens. 

• In particular, the results for automated delamination propagation analysis 
under quasi-static loading showed the following: 

o Using the default release tolerance (reltol= 0.2) as suggested in 
the Abaqus/Standard R handbook or increasing the value, as 
suggested in the user's manual, may help to overcome 
convergence problems, however, it leads to an undesired 
overshoot of the computed result compared to the benchmark, 
o A combination of release tolerance and contact stabilization is 
required to obtain more accurate results, 
o A gradual reduction of the release tolerance and contact 
stabilization over several analyses is suggested, 
o Good agreement between analysis results and the benchmarks 
could be achieved for release tolerance values (relto/< 0.1) in 
combination with contact stabilization (cs=lxl0~ 6 ). 

Overall, the results are promising and the current findings concur with previously 
published conclusions [3,4,6]. However, further assessment for mixed-mode 
delamination fatigue onset and growth is required. Additional studies should also 
include the assessment of the propagation capabilities in more complex mixed-mode 
specimens and on a structural level. 

Assessing the implementation in one particular finite element code illustrated the 
value of establishing benchmark solutions since each code requires specific input 
parameters unique to its implementation. Once the parameters have been identified, 
they may then be used as a starting point to model delamination growth for more 
complex configurations. 
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